% by ylx,clq on 6/15
% Version 2.0
%% 画图部分需要改进
%% for fig 10 12 14

function [prctages, cost_struct, profits_prct] = figure1x_even(j_all_all)
   
    global I J ci cI r CURRENT s_time
    ci_origin = ci;
    cI_origin = cI;
    r_origin = r;
   

    cI = 100;
    r = 100;
    cis = 40:10:90;
    N = size(j_all_all, 2); % 可否这样写省去N参数?
    len_seq = size(j_all_all, 1); %认为j_all_all为seq行向量的组合
    len_cis = size(cis, 2);
    index = zeros(1,len_seq); % 每个seq最优利润病人的总数
    cost_struct = zeros(3,I); %收入减去两个成本的矩阵
    opt_profits = zeros(1,N);  %同个seq的最高利润组合
    opt_costruct = zeros(3,N);  %所对应的cost structure
    max_profits = zeros(1,len_seq); %不同seq的最高利润
    max_costructs = zeros(3,len_seq); % 所对应的cost structure
    profits_prct = zeros(3,len_cis); % percentiles 对应的平均利润
    costructure = zeros(3,len_cis);
    i_star = zeros(1,N);
    slot_patientno = zeros(I,len_seq);
    prctages = zeros(I, len_cis);


    for jj = 1:len_cis
        % 这里大段代码重复 未来可以打包成函数
        % 代码来自 prop_policy.m
        ci = cis(jj);    
        for ii = 1:len_seq
            fprintf("\t fig1x_even -> %d/%d * %d/%d\n", jj, len_cis, ii, len_seq);
            if mod(ii, 10) == 0
                fprintf("\t\t Totle time = %.1f\n", ftimer()-s_time);
            end 
            Line = j_all_all(ii, :);
            S = zeros(I, J);
            for n = 1:N 
                CURRENT = n; % 用CURRENT控制了findQ的行为
                j = Line(n);
                profits = zeros(1, I);
                for i = 1:I
                    S_test = S;
                    S_test(i, j) = S_test(i, j) + 1; 
                    Q = findQ(S_test); % size(Q) = [I, CURRENT]
                    R = findR(Q);
                    [profit,~,~,P_Part] = function_fQR(Q, R); % P_part = parts of profits
                    profits(i) = profit;
                    %display(P_Part)
                    %return;
                    cost_struct(:,i) = P_Part(:,n);

                end
                i = find(profits == max(profits));
                i = i(1);
                S(i, j) = S(i, j) + 1;
                Q = findQ(S);
                i_star(n) = i;
                opt_profits(n) = profits(i);
                opt_costruct(:,n) = cost_struct(:,i);
                %display(opt_profits)
                % display(S)
                %return;
                % R = findR(Q); % R看起来对这个程序没用
            end
            % display(i_star);
            %return;
            max_profits(ii) = max(opt_profits); 

            index(ii) = find(opt_profits == max(opt_profits)); % 把index改为index(ii)
            max_costructs(:,ii) = opt_costruct(:,index(ii)); % 把index(1)改为 index(ii)
            for i = 1:I
                slot_patientno(i,ii) = sum(i_star(1:index(ii)) == i); % 每个slot分配的病人总数
                %disp(slot_patientno)
            end
            %percentages(:, jj) = percentages(:, jj) + sum(S, 2)/sum(sum(S)); % 先循环叠加再求占比和平均，可以减小误差
            
        end
        %display(slot_patientno)
        %display(index)
        % return;
        % display(max_profits)
        % display(prctile(max_profits,[10; 50; 90]))
        % return;
        prctages(:,jj) = sum(slot_patientno, 2)/sum(index);
        profits_prct(:, jj) = prctile(max_profits,[10; 50; 90]);
        costructure(:, jj) = mean(max_costructs,2);
        
    end
    %percentages = percentages / len_seq;

    % display(profits_prct)
    % display(costructure)
    % display(percentages)
    % display(cis)
    % >>> Need a Ploter Here
    figure(10); 
    for i = 1:I
        if mod(i,2)==0
            plot(cis, 100*prctages(i, :), '--',"LineWidth",2);
        else
            plot(cis, 100*prctages(i, :),"LineWidth",2);
        end
        hold on
    end
    xlabel("c_i");
    ylabel("Percentage of patients assigned");
    yticks([6:2:18]);
    xticks([40:10:90]);
    legend('1','2','3','4','5','6','7','8','Location','Best')
    ylim([6 18])
    % DRAWBRACE([88,11.5],[88,12.3],5) 网上查到可以，但是用起来不行
    %text(87,13,'\downarrow 2-7')
    %text(87,16,'\downarrow 1')
    %text(87,14,'\uparrow 8')
    grid on
    % <<< 

    % Fig 12
    figure(12);
    for i = 1:3
        plot(cis, profits_prct(i, :),"LineWidth",2);
        hold on
    end
    xlabel("c_i");
    ylabel("Expected Profit");
    xticks([40:10:90]);
    yticks([1000:50:1450]);
    legend("10 prctile","50 prctile" ,"90 prctile","location","Best");
    ylim([1000 1450]);
    grid on
    %input("Press ENTER to continue.\n");


    % Fig 14
    figure(14);
    for i = 1:3
        plot(cis,costructure(i,:),"LineWidth",2);
        hold on
    end
    xlabel("c_i");
    ylabel("Components of expected profit for various c_i");
    xticks([40:10:90]);
    yticks([0:500:2500]);
    legend("r E[\Sigma_iX_i]","c_i E[\Sigma_iY_i]" ,"c_I E[\Sigma_IY_I]","location","Best");
    ylim([0 2500]);
    grid on
    %input("Press ENTER to continue.\n")


    % >>> 恢复global变量的值（如果后续没用到原来的值，不写也行）
    ci = ci_origin;
    cI = cI_origin;
    r = r_origin;
    % <<<
    return;

